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Abstract. We consider the initial boundary value problem for the ho- 
mogeneous time-fractional diffusion equation dtu — Au — (0 < a < 
I) with initial condition u(x, 0) = v(x) and a homogeneous Dirich- 
let boundary condition in a bounded polygonal domain £1. We shall 
study two semidiscrete approximation schemes, i.e., Galerkin FEM and 
lumped mass Galerkin FEM, by using piecewise linear functions. We 
establish optimal with respect to the regularity of the solution error 
estimates, including the case of nonsmooth initial data, i.e., v £ La(fl). 



1. Introduction 

We consider the model initial-boundary value problem for the fractional 
order parabolic differential equation (FPDE) for u(x,t): 

dfu - Au = f(x, t), in ft T > t > 0, 

(1.1) u = 0, ondn T>t>0, 

u(0) = v, in ft 

where ft is a bounded polygonal domain in M d (d = 1, 2, 3) with a boundary 
dQ and v is a given function on £1 and T > is a fixed value. 

Here dfu (0 < a < 1) denotes the left-sided Caputo fractional derivative 
of order a with respect to t and it is defined by (see, e.g. [121 P- 91] or [221 
p. 78]) 

i r\. 



where T(-) is the Gamma function. Note that if the fractional order a tends 
to unity, the fractional derivative df u converges to the canonical first-order 
derivative ^ [12] , and thus the problem ( |1.1[ ) reproduces the standard par- 



abolic equation. The model (1.1) is known to capture well the dynamics of 
anomalous diffusion (also known as sub-diffusion) in which the mean square 
variance grows slower than that in a Gaussian process [I], and has found a 
number of important practical applications. For example, it was introduced 
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by Nigmatulin [21] to describe diffusion in media with fractal geometry. A 
comprehensive survey on fractional order differential equations arising in 
dynamical systems in control theory, electrical circuits with fractance, gen- 
eralized voltage divider, viscoelasticity, fractional-order multipoles in elec- 
tromagnetism, electrochemistry, and model of neurons in biology is provided 
in [3]; see also [52] . 

The capabilities of FPDEs to accurately model such processes have gen- 
erated considerable interest in deriving, analyzing and testing numerical 
methods for solving such problems. As a result, a number of numerical 
techniques were developed and their stability and convergence were inves- 
tigated, see e.g. O EU QU [201 [27] . Yuste and Acedo in [27] presented 
a numerical scheme by combining the forward time centered space method 
and the Grunwald-Letnikov method, and provided a von Neumann type sta- 
bility analysis. By exploiting the variational framework introduced by Ervin 
and Roop, |9j, Li and Xu [15] developed a spectral approximation method 
in both temporal and spatial variable, and established various a priori error 
estimates. Deng [6] analyzed the finite element method (FEM) for space- 
and time-fractional Fokker-Plank equation, and established a convergence 
rate of 0(T 2 ~ a + h 11 ), with a G (0, 1) and \i G (1, 2) being the temporal and 
spatial fractional order, respectively. 

In all these useful studies, the error analysis was done by assuming that 
the solution is sufficiently smooth. The optimality of the established es- 
timates with respect to the smoothness of the solution expressed through 
the problem data, i.e., the right hand side / and the initial data v, was not 
considered. Thus, these studies do not cover the interesting case of solutions 
with limited regularity due to low regularity of the data, a typical case for 
inverse problems related to this equation; see e.g., [3], [23] Problem (4.12)], 
and also (TO] [H] for its parabolic counterpart. 

There are a few papers considering construction and analysis of numerical 
methods with optimal with respect to the regularity of the solution error 
estimates for the following equation with a positive type memory terms 

1 /■* 

(1.2) d t u - —— / (t - r) a_1 Au(r)dr = f(x, t), t > 0, < a < 1, 

r(a) Jo 



This equation is closely related, but different from (1.1). For example, 
McLean and Thomee in \17\ [TB] developed a numerical method based on 
spatial finite element discretization and Laplace transformation with quadra- 



tures in time for (1.2) with a homogeneous Dirichlet boundary data. In [17} 
Theorem 5.1] the convergence of the proposed method has been studied and 
maximum-norm error estimates of order 0(i _1_a /i 2 £|), £h = |l n ^-|) were 
established for initial data v G Lo^fi). Further, in [18^ Theorem 4.2] a 
maximum-norm error estimate of order 0{h?(^) was shown for smooth ini- 
tial data v G H 2 . Mustafa [20] studied a semidiscrete in time and and fully 
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discrete schemes, Crank-Nicolson in time and finite elements in space, and 
derived error bounds for smooth initial data; see, e.g. [20^ Theorem 4.3]. 

The lack of optimal with respect to the regularity error estimates for the 
numerical schemes for FPDEs with nonsmooth data is in sharp contrast with 
the finite element method (FEM) for standard parabolic problems, a = 1. 
Here the error analysis is complete and various optimal with respect to the 
regularity of the solution estimates are available |25| . The key inngredient 
of the analysis is the smoothing property of the parabolic operator and its 
discrete counterpart |25| Lemmas 3.2 and 2.5]. For the FPDE (1.1), such 
property has been established recently by Sakamoto and Yamamoto |23j ; 
see Theorem 12.11 below for details. 

The goal of this note is to develop an error analysis with optimal with 
respect to the regularity of the initial data estimates for the semidiscrete 
Galerkin and the lumped mass Galerkin FEMs for the problem (1.1) on 
convex polygonal domains. 

Now we describe our main results. We shall use the standard notations 
in the finite element method [25] . Let {7/i} 0<ft<1 be a family of regular 
partitions of the domain f2 into d-simplexes, called finite elements, with h 
denoting the maximum diameter. Throughout, we assume that the triangu- 
lation Th is quasi-uniform, that is the diameter of the inscribed disk in the 
finite element r G Th is bounded from below by h, uniformly on Th- The ap- 
proximate solution Uh will be sought in the finite element space Xh = Xft(J2) 
of continuous piecewise linear functions over the triangulation Th 

X is a linear function over r, Vr G Th} ■ 

The semidiscrete Galerkin FEM for the problem (1.1) is: find Uh(t) G Xh 
such that 

(d?u h ,x) + a(u h ,x) = (J,x), VxE4T>(>0, 

Uh{0) = V h , 

where a(u,w) = (Vu,Vw) for u, w G Hq(Q), and Vh G Xh is a given 
approximation of the initial data v. The choice of Vh will depend on the 
smoothness of the initial data v. Following Thomee [25], we shall take 
Vh = RhV in case of smooth initial data and Vh = PhV in case of nonsmooth 
initial data, where Rh and Ph are Ritz and the orthogonal L2(f2)-projection 
on the finite element space Xh, respectively. 

We shall study the convergence of the semidiscrete Galerkin FEM (1.3) 
for the case of initial data v G H q (Q), q = 0,1,2 (for the definition of these 
spaces, see Section 2.2). The case q = 2 is referred to as smooth initial data, 
while the case q = is known as nonsmooth initial data. 

In the past, the initial value problem for a standard parabolic equation, i.e. 
a = 1, has been thoroughly studied in all these cases. It is well known that, 
for smooth initial data, the solution Uh satisfies an error bound uniformly 
in t > [25, Theorem 3.1]: 

(1.4) || uft (t) -u(t)\\ +h\\V(u h (t) -u(t))\\ < Ch 2 \\v\\ 2 , for t > 0. 
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We also have a nonsmooth data error estimate, for v assumed to be only in 
£2(0), but which deteriorates for t approaching |25} Theorem 3.2], namely, 

(1.5) \\u h (t) -u(t)\\ +h\\V(u h (t) -u(t))\\ < Ch 2 ^ 1 ^, for t > 0. 

The proof of all these results exploits the smoothing property of the parabolic 
problem via its representation through the solution operator E(t) = e tz \ 
namely, 

u(t) = E(t)v+[ E(t- s)f(s)ds, i>0. 
J 

In this paper we establish analogous results for the semidiscrete Galerkin 



FEM (1.3) for the model problem (1.1). The main difficulty in the error 



analysis stems from limited smoothing properties of the FPDE, cf. Theo- 



rem 



2.1 Note that the solution operator for the FPDE is defined through 



the Mittag-Leffler function, which decays only linearly at infinity, cf. Lemma 



2.1, in contrast with the standard parabolic equation whose solution decays 
exponentially for t — > 00. The difficulty is overcome by exploiting the map- 
ping property of the discrete solution operators. 

Our main results follows. Firstly, in case of smooth initial data, we 



derived the same error bound (1.4) uniformly in t > (cf. Theorem 3.1), 



as is in the case of the standard parabolic problem. Secondly, for quasi- 
uniform meshes we derived a nonsmooth data error estimate, for v £ Li2(P) 



only, which deteriorates for t approaching (cf. Theorem 3.2) 

(1.6) \\u h {t)-u(t)\\+h\\X7(u h (t)-u{t))\\ < Ch 2 £ h t- a \\v\\, £ h = \hxh\,t> 0. 

This result is similar to the counterpart of standard parabolic problem but 
derived for quasi- uniform meshes and with an additional log- factor, i^. 

Further, we study the more practical lumped mass semidiscrete Galerkin 
FEM. We have shown the same rate of convergence for the case of smooth 



initial data (cf. Theorem 4.1), and also the almost optimal error estimate 
for the gradient in the case of data v G H l {^l) and v G L^{^1) (see estimate 



(4.2)). For the case of nonsmooth data, v E £2(0), f° r general quasi-uniform 



meshes, we were only able to establish a suboptimal L2-error bound of order 
0(h£ht~ a ), see (4.15). Further, inspired by the study in [2], we also consider 



special meshes. Namely, for a class of special triangulations satisfying the 



condition (4.16), which holds for meshes that are symmetric with respect to 



each internal vertex [21 Section 5], we show an almost optimal convergence 
estimate (4.17): 

\\u h (t)-u(t)\\ <Ch 2 £ h t- a \\v\\, 

where Uh(t) is the solution of the lumped mass FEM. This estimate is similar 
to the one derived for the lumped mass semidiscrete Galerkin method for 
the standard parabolic equation [21 Theorem 4.1]. 



Finally, in Theorem 5.1 we establish a superconvergence result for the 



postprocessed gradient of the error in case of smooth initial data and a 
planar domain for special meshes. This improves the convergence order in 
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.ff^-norm from 0(h) to 0(h 2 lht~?) for both Galerkin and the lumped mass 
finite element approximation. 

The paper is organized as follows. In Section [2l we state basic properties 



of the Mittag-Leffler function, the smoothing property of the equation ( 1.1 ), 
and some basic estimates for finite element projection operators. In Sections 
[3] and [4j we derive error estimates for the standard Galerkin FEM and 
lumped mass FEM, respectively. In Section [5] we give a superconvergence 
result for the gradient of the error in case of smooth initial data. Finally, in 
Section [6] we present some numerical tests for a number of one-dimensional 
examples, including both smooth and non-smooth data. The numerical tests 
confirm our theoretical study. 

We assume that the mesh size h of the triangulation Th satisfies < h < 1. 
Throughout we shall denote by C a generic constant, which may differ at 
different occurrences, but is always independent of the mesh size h, the 
solution u and the initial data v. 

2. Preliminaries 
In this section, we collect useful facts on the Mittag-Leffler function, reg- 



ularity results for the fractional diffusion equation ( 1.1 ), and basic estimates 
for the finite-element projection operators. 

2.1. Mittag-Leffler function. We shall use extensively the Mittag-Leffler 
function E a ^(z) defined below 

00 y. 

where T(-) is the standard Gamma function defined as 

/•oo 

r(z) = / f- 1 e- t dt $t(z) > 0. 
J 

The Mittag-Leffler function is a two-parameter family of entire functions 
of z of order a -1 and type 1 [121 pp. 42]. The exponential function is 
a particular case of the Mittag-Leffler function, namely E\ t \(z) = e z , \12\ 
pp. 42]. Two most important members of this family are E a i(— Xt a ) and 
* a ~ 1 -E«,a(-At a ), which occur in the solution operators for the initial value 



problem and the nonhomogeneous problem (1.1), respectively. There are 
several important properties of the Mittag-Leffler function E a ^p(z), mostly 
derived by M. Djrbashian (cf. [7, Chapter 1]). 

Lemma 2.1. Let < a < 2 and (3 G M be arbitrary, and ^ < \i < 
min(7r,a7r). Then there exists a constant C = C(a,f3,fi) > such that 

(2.1) \E*A*)\ ^ T7TT ^ - l arg ^l - w - 

~r I % I 

Moreover, for A > 7 a > 0, and t > we have 

(2.2) dfE a) i(-Xt a ) = -\E aA (-\t a ). 
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Proof. The estimate (2.1) can be found in [12], pp. 43, equation (1.8.28)] 



or |22] Theorem 1.4], while (2.2) is discussed in |12] Lemma 2.33, equation 



(2.4.58)]. □ 
2.2. Solution representation. To discuss the regularity of the solution of 



(1.1), we shall need some notation. For q > 0, we denote by H q (Q) C L%{£l) 



the Hilbert space induced by the norm 



with (•, •) denoting the inner product in Li2(P) and {\j}JL l and { t fj}JL 1 be- 
ing respectively the eigenvalues and eigenfunctions of —A with homogeneous 
Dirichlet boundary data on d£l. The set {<^j}^ =1 forms an orthonormal ba- 
sis in £2(0). Thus \v\q = \\v\\ = {v^v) 1 / 2 is the norm in L2(Q), \v\i the 
norm in Hq = Hq(Q) and | t; 1 2 = ||Au[| is equivalent to the norm in H 2 (£l) 
when v = on dQ, [25]. We set H~ q = (H q )' ', the set of all bounded linear 
functionals on the space H q . 



Now we give a representation of the solution of problem (1.1) using the 



Dirichlet eigenpairs {(Xj, fj)}- First, we introduce the operator E(t): 

00 

(2.3) E(t)v = Y j E a>1 {-\ j t a ) {v )ipj )ipj{x). 



3=1 



This is the solution operator to problem ( |1.1[ ) with a homogeneous right 
hand side, so that for f(x, t) = we have u(t) = E(t)v. This representation 



follows from an eigenfunction expansion and (2.2) [23] . Further, for the non- 
homogeneous equation with a homogeneous initial data v = 0, we shall use 
the operator defined for \ G L 2 (Q) as 



(2.4) E(t) X = 2^t a - l E aia (-\ j t a ) (x,<Pj) <Pj(x). 

j=0 

The operators E(t) and E{t) are used to represent the solution u(x, t) of 



(]Tlj): 

t 



u(x, t) = E(t)v + I E(t- s)f(s)ds. 
Jo 



It was shown in [23] Theorem 2.2] that if f(x,t) £ L 2 ((0, T); L 2 (0)), then 
there is a unique solution u(x,t) € ^((0, T); H 2 (£l)). 



For the solution of the homogeneous equation ( 1.1 ), which is the object of 
our study, we have the following stability and smoothing estimates, essen- 
tially established in [23] Theorem 2.1], and slightly extended in the theorem 
below. Since these estimates will play a key role in the error analysis of the 
FEM approximations, we give some simple hints of the proof. 
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Theorem 2.1. The solution u(t) 
satisfies the following estimates 



E(t)v to problem (1.1) with / = 



(2.5) 



\{d?) l u(t)\ p <Cr^ l+ ^\v\ q , t>0, 



where for I = 0, < g < p < 2 and for 1 = 1, < p < q < 2 and q < p + 2. 

Proof. First we discuss the case i = 0. According to parts (i) and (iii) of 
Theorem 2.1], we have 

< cr a(1 -^\v\ 



(2.6) 



\u{t)\ 2 + \\d?u 



19) 



q = 0,2. 



By means of interpolation of estimates (2.6) for q = and q = 2, we get the 



desired estimate (2.5) for the the case p = 2, < g < 2. 



Further, applying part (i) of |23t Theorem 2.1], we have 
(2.7) \\u(t)\\ <C\\v\\. 



Thus, interpolation of (2.6) for q = 2 and (2.7) yields (2.5) for < p = q < 2. 



The remaining cases, < q < p < 2, follow from the interpolation of (2.6) 
with q = and (2.7). This shows the assertion for i = 0. 

1. It follows from the representation formula 



Now we consider the case 



(2.3) and Lemma |2.1| that 

oo 



\d?u{t)\l 



oo 

= r«( 2+ ^ ^(A,t a ) 2+p -^ aj i(-A^ a ) 2 Aj(«,^; 



sup 



^5 (1 + A^°) 2 



^A 9 (w,99j) 2 < Ct~ a( - 2+p ~ q ^\ 



where we have used the fact that, in view of Young's inequality, 
{\ j t a ) 2+ P- ( i 



3 (1 + XjV*)* 



<C for p<q<p + 2. 



Thus, we get 



\d?u(t)\ p < Ct- a(1+ V)| 



This completes the proof of the theorem. 



□ 



Remark 2.1. Note that for I = 1 we have the restriction p < q, which is 
not present in the similar result for the standard parabolic problem, see, e.g. 
|25} Lemma 3.2]. This reflects the fact that FPDE has limited smoothing 
properties. The limited smoothing is also valid for the semidiscrete Galerkin 



approximation (see Lemma 3.1), which will influence the error estimates for 
the finite element solution. 
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2.3. Properties of Ritz and L2-projections on Xh- In our analysis we 
shall also use the orthogonal ^-projection Ph : ^(O) —> Xh and the Ritz 
projection : Hq(Q) — > Xh defined by 

respectively. It is well-known that the operators Ph and Rh have the follow- 
ing approximation properties. 

Lemma 2.2. The operators Ph and Rh satisfy 

(2.9) \\P h ^-^\\+h\\V{P h ^-^)\\ <Ch*\i/>\ q , for i>eH«, q = 1,2. 

(2.10) 11^-^11+ h\\V(R h rP-^)\\ <Ch q \iP\ q , for ^€H q ,q = l,2. 
In particular, (|2.9l) indicates that Ph is stable in H 1 . 



Proof. The estimates (2.10) are well known, cf. e.g. |25|, Lemma 1.1] or 
[H Theorem 3.16 and Theorem 3.18]. For globally uniform meshes, the case 
considered in this paper, the H 1 stability of Ph directly follows from the error 



bound (2.9) and the inverse inequality. However, for more general meshes 
such stability is valid only under some mild assumptions on the mesh; see, 
e.g. @]. □ 

3. Semidiscrete Galerkin FEM 

In this section we derive error estimates for the standard semidiscrete 
Galerkin FEM. First we recall some basic known facts for the spatially 
semidiscrete standard Galerkin FEM. We begin with the smoothing prop- 
erties of the solution operators for the semidiscrete method as well as other 
preliminary results needed in the sequel. The error estimates hinge crucially 



on the smoothing properties of the discrete operator Eh, cf. (3.3). 



3.1. Semidiscrete Galerkin FEM and its properties. Upon introduc- 
ing the discrete Laplacian A/, : Xh — > Xh defined by 

-(A h v, x) = W, Vx) W>, x e X h , 

and fh = Phf we may write the spatially discrete problem ( |1.3| ) as 
(3.1) dfu h {t) - A h u h (t) = f h (t) for t>0 with u h (0) = v h . 



Now we give a representation of the solution of (3.1) using the eigenvalues 



and eigenfunctions {Xj}j =1 and {(fj}j =t of the discrete Laplacian —A 



h- 



First we introduce the discrete analogues of (2.3) and (2.4) for t > 0: 

N 

(3.2) E h (t)v h = Y J E a ,x{-X h j t a ){v h , 

3=1 
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N 

(3.3) E h (t)f h =J2t° l ~ 1 Ea,a(-fy a ) (fh,<Pj)<fy 

3=1 



Then the solution Uh(x, t) of the discrete problem (3.1 ) can be expressed by: 

(3.4) u h (x,t) = E h (t)v h + [ E h (t- s)f h (s)ds. 

Jo 

Also, on the finite element space X^, we introduce the discrete norm ||| • 
for any p £ R defined by 

N 

(3-5) IIMIi; = £( A *W'¥#) a ^4 

i=i 

Since we are dealing with finite dimensional spaces, the above norm is well 
defined for all real p. From the very definition of the discrete Laplacian — 
we have IHV'llli = IV'li an d also IHV'lllo = IIV'II f° r anv ^ £ -^Ti- So there will 
be no confusion in using \ijj\ p instead of IHV'lllp for p = 0, 1 and ip £ X^. 

To analyze the convergence of the semidiscrete Galerkin method, we shall 
need various smoothing properties of the operator E^t), which are discrete 
analogues of those formulated in (2.5). The estimates will be used for ana- 
lyzing the convergence of the lumped mass FEM in Section |4| 

Lemma 3.1. Let Eh(t) be defined by ( |3.2[ ) and 6 Xh- Then 

(3.6) \\\(d?y Uh (t)\\\p = \\\m £ E h (t)vh\\\ P < Ct-^ + ^\\\v h \\\ q , t > 0, 
where for I = 0, q < p and < p — q < 2 and for £=l,p<q<p + 2. 

Proof. First, consider the case £ = 0. Then using the representation ( |3.2[ ) of 
the solution Uh(t) and the bound for the Mittag-Leffler function E a a{z) in 
Lemma 2.1 we get for q < p 

N N 

\\\Mt)\\\l = 5>^iMt),^)i 2 = ^j) p \E a ,i(-^n\ 2 \(v h ,^)\ 2 

3=1 3=1 

N (X h t a ) p ~ q 



N 

<Ct- a ^«) J^(Aj)»|(« ft ,^)| 2 = C't- 0, ^)||| V/l |||2. 

3=1 



Here in the last inequality we have used the fact that for q < p and p < q < 
p + 2 the expression m&Xj(\jt a ) p ~ q / (1 + X!-t a ) 2 is bounded. 
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The estimates for the case I = 1 are obtained analogously using the 



representation (3.2) of the solution Uh{t) for p < q < p + 2 and Lemma 2.1 

N 

\\\dr Uh (t)\\\i = Y,(^r\(oru h (t),^)\ 2 

3=1 
N 

= j2(^) 2+ nE a A-\^n\ 2 \(v h ,^)\ 2 . 

5=1 



Now using the bound of Mittag-Leffler function in Lemma 2.1 and Young's 
inequality, we obtain 

= c , t -(2 Q! +a ! (p- 3 ))|||^|||2_ 

The desired estimate follows from this immediately. □ 

The following estimates are crucial for the a priori error analysis in the 
sequel. 



Lemma 3.2. Let E^ be defined by (3.3) and ip 6 X^. Then we have for all 
t > 0, 



(3.7) \\\E h (t)i;\\\ p < 



Ct-^+^MWg, p-2<q<p, 
Ct- 1+a M\\ q , p<q. 



Proof. By the definition of the operator Eh(t) and using Lemma 2.1 for 
E a ,a(z), we have for any p£l and q < p 

N 

(\ h t a \p-i N 

< ct -2 +a{2+q - P) mfX <_^_> £ (Ajm J]? 

\ j I j=l 

where for getting the last inequality we took into account < p — q < 2. 
The desired assertion for p < q follows from the fact that the eigenvalues 
{A^} are bounded away from zero independently of the mesh size h. □ 

Remark 3.1. Lemma \3.S\ expresses the smoothing properties of the operator 
Eh. While p = 0, 1, the parameter q can be arbitrary as long as it complies 
with the condition p— 2 < q < p. This flexibility in the choice of q is essential 
for deriving error estimates for problems with initial data of low regularity. 
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Further, we shall need the following inverse inequality. 

Lemma 3.3. There exists a constant C independent of h such that for all 
ip £ Xh we have for any real I > s 

\i < ch s - l \ 



(3.8) 

Proof. For quasi-uniform triangulations 7~h the inverse inequality \ip\i < 
Ch~ 1 \\ip\\ holds for all ?p G X^. By the definition of — this implies 
maxi<j<AT < C/h 2 . Thus, for the norm ||| • ||| p defined in (3.5), we obtain 
that for any real / > s 



N 



° i=i 
That completes the proof. 



□ 



3.2. Error estimates for smooth initial data. Here we establish error 
estimates for the semidiscrete Galerkin method for initial data v £ H 2 (Q). 
In a customary way we split the error Uh{t) — u(t) into two terms as 



(3.9) 



u h -u = (u h - R h u) + (R h u - u) := •& + g. 



By (2.10) and (2.5) we have for any t > and q = 1,2, 



(3.10) 



\\ g (t)\\ + h\\Ve(t)\\<ChH 



so it suffices to get proper estimates for $(t), which is done in the following 
lemma. 



Lemma 3.4. Let u and Uh be the solutions of (1.1) and (1.3), respectively, 
with Vh = RhV. Then for $(t) = Uh(t) — Rhu{t) we have 

(3.11) ||tf(t)|| + h||W(i)|| < Ch 2 \v\ 2 . 
Proof. We note that i? satisfies 

(3.12) fl«0(i) - A h t?(t) = -P h d?g{t) for t > 0. 

For v G H q , q = 1, 2 the Ritz projection R^v is well defined, so that i?(0) = 
and hence, by Duhamel's principle (3.4), 



(3.13) 



E h (t-s)P h d? g (s)ds. 







By using Lemma 3.2 with p = 1 and q = 0, the stability of Ph, (2.10), and 
the estimate (2.5) with £ = 1, p = 1 and we find, for q = 1, 2, 

(3.14) 



VE h (t - s)P h d? Q (s)\\ < C{t - s)" 1 Wq{s)\\ 



< Ch(t - s)^' 1 Idfuis)^ 

< Ch{t-s)^- x s a ^ + ^\v 
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By substituting this inequality into (3.13) we obtain that for q = 1, 2 
(3.15) ||V0(t)|| <Ch (t- s )f" V(-3+§) ds \v\ q < Chr a ^-^\v\ q , 



o 



where we have used that for a < 1 



(i_ s) f-V(-!+!) ds 



a Set I ga 

^ 2 2 2 



(1 



S )f-i s «(-|+f)d s 



fld.aC-l + fJ + l)^ 1 -*) 



with B(-, •) being the standard Beta function. Since both arguments, f > 
and z ^ 2 ±1 a + 1 > for q = 1, 2, the value B(a,a(-§ + |) + 1) is finite. 
Taking q = 2 we get the desired estimate for Vi9. 



Next, by using the smoothing property of the operator Eh in Lemma 3.2 
with p = q = and that of the operator E in Theorem 2.1 with 1 = 1 and 
p = q = 2, we get 

< / \\E h {t-5)P h d? e (s)\\ds 



(3.16) 



o 

< C I (t-s 

Jo 

< Ch 2 [ (t- 

Jo 

< Ch 2 f (t- 

Jo 



i«- 1 [|S? e ( a )[|d a 

„\a— 1 „— a 



s Q ds|t>|2 = CB(a, 1 — a)/i |v (2- 



This completes the proof. 



□ 



Using the triangle inequality and the estimates (3.10) and (3.11) we get 
the main result in the subsection. 



Theorem 3.1. Letu anduh be the solutions of (1.1) and (1.3), respectively, 
with Vh = RhV- Then 

(3.17) \\u h (t) - u(t)\\ + h\\X7(u h (t) - u(t))\\ < Ch 2 \v\ 2 . 



Remark 3.2. As a byproduct of estimates (2.10) and (3.15) we also got a 



bound for the error for v G H (Q) and Vh = Rh v: 

(3.18) l|VK(i) - u(i))|| < Cht~^\v\i. 

Remark 3.3. I n vi ew of the smoothing property of the operator Eh estab- 
lished in Lemma 3.2, we can improve the estimate ofd(t) for q = 2 to 0(h 2 ) 
at the expense of slightly increasing the factor by 0(t~z): 

\\VE h {t - s)P h d?Q{s)\\ < Ch\t - s)!- 1 \d?u(s)\ 2 

< Ch 2 (t- s)f-V>| 2 , 

which yields 

(3.19) ||W|| < CTi 2 * - ? \v\ 2 . 
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3.3. Error estimates for non-smooth initial data. Now we prove an 
error estimate for nonsmooth initial data, v S £2(0), an d the intermediate 
case, v £ Since the Ritz projection R^v is not defined for v 6 L2(£l), 

we shall use instead the ^-projection Ph onto the finite element space Xh, 
and split the error Uh — u into: 

u h -u = (uh- Phu) + {PhU - u) := ■& + q. 



By Lemma 2.2 and Theorem 2.1 we have 

(3.20) ||e(t)|| + h\\Vg{t)\\ < Ch 2 \u(t)\ 2 < Ch 2 r a ^~2 ) \\v\\ q , q = 0,1. 

Thus, we only need to estimate the term Obviously, Phdfg = dfPh(PhU- 
u) = and we get the following problem for 

(3.21) d?$(t)-A h d(t) = -A h (R h u-P h u)(t), t>0, 0(0) = 0. 



Then with the help of formula (3.3), $(t) can be represented by 



(3.22) 



0(t) = - / E h (t - s)A h (R h u - P h u)(s) ds. 
Jo 



Next, we show the following estimate for -d(t): 



Lemma 3.5. Let •&{€) be defined by (3.22). Then for p = 0, 1, q = 0, 1, and 



lh = I hxh\, the following estimate holds 

\W(t)\\ P <ch 2 -n h t- a ^\\v\\ q . 

Proof. By Lemma 3.2 with p = 0, 1 and q = p — 2 + e, for any e > we have 

\\m\\p < f \\E h (t - s)A h (R h u - P h u)(s)\\ p ds 
Jo 

< [ (t- S )i Q - 1 |||A^( J R^-P h n)||| p _ 2+e d S 
Jo 

< [ (t - s)^ a ~ 1 \\\R h u - P h u\\\ p+e ds := A. 
Jo 



Further, we apply the inverse inequality (3.8) for R^u — -f\u, the bounds 
(2.9) and (2.10), for P^u — u and R^u — u, respectively, and the smoothing 
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property (2.5) with £ 



A < ChT 



and p = 

t 

(t-a)3 



2 to get 
a ~ l \\R h u 



o 



P h u\\ p ds 

u(s)\\2<ls 



< Ch 2 - p ~ e / (t 
Jo 
ft 

Jo 



CB{-a,l 
C 



a H — a ) h 



,2-p- ei -a(l 



-f) 



< — fc 2 -P- e t-«( 1 -t) 



It/- 



The last inequality follows from the fact -B(|a, 1 — a + |a) — — ^g+e^ 

and T(|a) ~ ~ as e — > + , e.g., by means of Laurenz expansion of the 
Gamma function. The desired assertion follows by choosing e = 1/lh- D 



Then Lemma 3.5 and the triangle inequality yield the following almost 



optimal error estimate for the semidiscrete Galerkin method for initial data 
v G H q , q = 0,l: 



Theorem 3.2. Let u and Uh be the solutions of (1.1) and (1.3) with Vh 
Ph,v, respectively. Then with 1^ = | ln/i| 

\V(u h (t)-u(t))\\ <Ch 2 £ h t- a ^-^ 



(3.23) \\u h (t)- 



+ i 



Remark 3.4. For v 6 H l (Q) and Vh 



q, 9 = 0, 1. 

RhV, we have established the esti- 



mate (3.18), which is slightly better than (3.23), since it does not have the 



factor 1^. 

3.4. Problems with more general elliptic operators. The preceding 
analysis could be straightforwardly extended to problems with more gen- 
eral boundary conditions/spatially varying coefficients. In fact this is the 
strength of the finite element method and the advantages of the direct nu- 
merical methods for treating such problems in comparison with some an- 
alytical techniques that are limited to constant coefficients and canonical 



domains. More precisely, we can study problem (1.3) with a bilinear form 
a(-, •) : V x V i-)- R of the form: 



a(u,x) = / (k(x)Vu ■ Vx + c(x)ux) dx, 



(3.24) 



where k(x) is a symmetric d x d matrix-valued measurable function on the 
domain Q with smooth entries and c(x) is an Loo-function. We assume that 

c |£| 2 < fk{x)i < ci|£| 2 , for (el 11 , iGfl, 

where cq,ci > are constants, and the bilinear form a(-,-) is coercive on 
V = iL 1 (Sl). Further, we assume that the problem a(w, x) = (f, x)> G V 



2.2 
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has unique solution w G V, which for / G £2^) has full elliptic regularity, 

\H\ H 2<c\\f\\ L2 . 

Under these conditions we can define a positive definite operator A : 
Hq — > H , which has a complete set of eigenfunctions <fj(x) and respective 
eigenvalues Xj(A) > 0. Then we can define the spaces H q as in Section 
Further, we define the discrete operator Ah ■ Xh — > Xh by 

{Ahi>, x) = aO, x), W>, x G V h . 

Then all results for problem ( |f can be easily extended to fractional-order 
problems with elliptic equations of this more general form. 

4. Lumped mass finite element method 

Now we consider the lumped mass FEM in planar domains (see, e.g. 
Chapter 15, pp. 239-244]). For completeness we shall introduce this 



approximation. Let zj, j = 1, . . . ,d + 1 be the vertices of the d-simplex 
t G Th- Consider the quadrature formula 

d+i 

(4.1) Q T , h (f) = j^y^ /(4) ~ / f dx - 



3=1 jt 



We may then define an approximation of the L2-inner product in Xh by 

(4.2) (w,x)h=^2Qr,h(w X )- 

Then lumped mass Galerkin FEM is: find Uh(t) G Xh such that 

s {d?u h ,x)h + a(u h ,x) = (f,x) VxGX h , t > 0, 

(4.3) , . 

V ' «fc(0) = v h . 

We now introduce the discrete Laplacian — A^ : Xh — > Xh, corresponding 
to the inner product (•, -)h, by 

(4.4) -(A/^,x)/» = W,Vx) W, X €X h . 
Also, we introduce the .^-projection, Ph : £2(^2) — ^ Xh by 

(P h f,x)h = (f,x), v x ei, 

The lumped mass method can then be written with fh = Phf in operator 
form as 

dfu h {t) - A h u h (t) = f h {t) for t > with u A (0) = v h . 
Similarly as in Section [3j we define the discrete operator Fh by 

N 

(4.5) F h (t)v h = E atl {-\*}t a )(v h , <p)) h q>), 

3=1 

where {Aj}jLi and x are respectively the eigenvalues and the or- 

thonormal eigenfunctions of — A^ with respect to (•,•)/»■ 
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Analogously to (3.3), we introduce the operator Fh by 

N 

(4.6) F h f h (t) = ^t a ~ 1 E a)a {-\ h j t a )(f h ,(p h j ) h (p h j . 

i=i 



Then the solution Uh to problem (4.3) can be represented as follows 



u h (t) = F h (t)v h + / F h (t - s)f h (s)ds. 
Jo 

For our analysis we shall need the following modification of the discrete 



norm (3.5), 
(4.7) 



| p , on the space Xh 
N 



The following norm equivalence result is useful. 



Lemma 4.1. The norm 



defined in (4.7) is equivalent to the norm 



| • \p on the space Xh for p = 0, 1. 

Proof. The equivalence the the two norms for the case of p = is well known: 

From the definitions of the discrete Laplacian — Ah and the eigenpairs {(A^ , tpj 
we deduce 



< 



< 



0) 



W> G X h . 



N 



HWH 2 = {{-A h )^i>) h = £aJ(v,#)£ 

i=i 

This completes the proof of the lemma. 



□ 



We shall also need the following inverse inequality, whose proof is identical 
with that of Lemma 13.31 



(4.8) 



< Ch 



s-l 



I > S. 



We show the following analogue of Lemma 3.2 



Lemma 4.2. Let Fh be defined by (4.6). Then we have for ip e Xh and all 

t > 0, 



\mmwp < 



Ct~ 1+a 



q i 



q , P~2<q<p, 
p < q. 



Proof. The proof essentially follows the steps of the proof of Lemma |3.2| by 
replacing the eigenpairs (A^ , tpj) by (Xj , (pj), and the L2-hiner product (•, •) 
by the approximate L2-inner product (-, -)h and thus it is omitted. □ 
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We need the quadrature error operator Qh : Xh — > Xh defined by 
(4.9) (VQ hX ,V^) =e h ( X ,1>) ■= M)h~M) V X ,V G X h . 

The operator Qh, introduced in [2], represents the quadrature error (due to 
mass lumping) in a special way. It satisfies the following error estimate: 



Lemma 4.3. Let Ah and Qh be the operators defined by (4.4) and (4.9), 
respectively. Then 

\\VQ h x\\ + h\\A h Q hX \\ < C7^ +1 ||VP X || V X £4p = 0, 1. 

Proof. See [21 Lemma 2.4]. □ 

4.1. Error estimate for smooth initial data. We shall now establish 
error estimates for the lumped mass FEM for smooth initial data, i.e., v E 
H 2 (Q). 



Theorem 4.1. Letu anduh be the solutions of (1.1) and (4.3), respectively, 
with Vh = RhV- Then 



\u h (t) - u(t)\\ + h\\V{u h {t) - u(t))\\ < Ch 2 \v 



2- 



Proof. Now we split the error into Uh(t) — u{t) = Uh(t) — u(t) + S(t) with 
S(t) = Uh(t) — Uh(t) and Uh(t) being the solution by the standard Galerkin 
FEM. Since we have already established the estimate ( 3.17| ) for Uh — u, it 
suffices to get the following estimate for S(t): 

(4.10) \\5(t)\\+h\\V5(t)\\ < Ch 2 \v\ 2 . 

It follows from the definitions of the Uh(t), Uh{t), and Qh that 

d?5{t) - A h S(t) = h h Q h dfu h (t) for t > 0, 5(0) = 

and by Duhamel's principle we have 

S(t) = f F h (t - s)A h Q h d?u h (s)ds. 
Jo 

Using Lemmas |4.1[|4.2[ and |4.3| we get for % 6 Xh' 

\\X7F h (t)A h Q hX \\ < Ct^WAhQhxW < Ct%- l h\\V X \\- 
Similarly, for x G Xh 

\\F h (t)A h Q hX \\ < Ctf^lHAftQhxIII-i < Ctf^WVQhxW < Ct^h^xl 
Consequently, using Lemma |3.1| with I = 1 , p = 1 and q = 2 we get 

||tf(t)|| + h\\V5(t)\\ < Ch 2 f (t - s^W^Uhis^ds 



o 

< Ch 2 ( (t-s)f- 1 s-t c / s |||^(0)||| 2 . 
Jo 

Since AhRh = PhA, we deduce 

IK(0)||| 2 = ||A ft iJ h «(0)|| = ||i^A«(0)|| < |«(0)| 2 < C|M| 2 , 



which yields (4.10) and concludes the proof. □ 
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An improved bound for ||V£(i)|| can be obtained as follows. In view of 
Lemmas 4.1 and 4.3 and (4.8), we observe that for any e > and x £ Xh 

\\VF h (t)A h Q hX \\ < C^ a - 1 |||A ft g & xl||-i+ e < Cif^^UVxII. 



Consequently, 

(4.11) ||V5(t)|| < Ch 2 ~ e 



(t 



\d?u h (a)\\\ida. 



o 



Now, to (4.11 ) we apply Lemma 3.1 with I = 1, p = 1 and q = 2 to get 

r-t 



\\V5(t)\\ < Ch 



2-e 



.1 







2 < C-h z ~H 
e 



Remark 4.1. In t/ie akwe estimate, by choosing e = 1/ih, 
get 

(4.12) ||W(t)||<C7i 2 4t" f M2, 



2 hb- 

| lnh\, we 



which improves the bound of ||V<5(t)|| for fixed t > by almost one order. 



Remark 4.2. Instead, if we apply to (4.11) Lemma 3.1 with I = 1, p = 1 
and q = 1 we shall get an improved estimate for 5(t) in the case of initial 
data v £ H 1 : 



(4.13) 



||V5(i)|| < Ch 2 e h t~ a \v\ 



4.2. Error estimates for nonsmooth initial data. Now we consider the 
case of nonsmooth initial data, i.e., v G £2^) as well as the intermediate 
case v G H l . Due to the lower regularity, we take = PhV. As before, 
the idea is to split the error into — u(t) = Uh(t) — u(t) + 6(t) with 

5(t) = Uh(t) — Uh(t) and Uh(t) being the solution of the standard Galerkin 
FEM. Thus, in view of estimate (3.23) it suffices to establish proper bound 
for S(t). 



Theorem 4.2. Let u andu^ be the solutions of (1.1) and (4.3), respectively, 
with Vh = PhV- Then with £h = \ lnh\, the following estimates are valid for 
t > 0: 



(4.14) 

and 

(4.15) 



\V(u h (t) - u(t))\\ < CM h t~ a{l -^\v\ q q = 0, 1, 



\u h (t) - u(t)\\ < CW+Hht-^-iMq q = 0,1. 



Furthermore, if the quadrature error operator defined by (4.9) satisfies 
(4.16) \\QhX\\ < Ch 2 \\ X \\ Vx £ Xh, 

then the following almost optimal error estimate is valid: 



(4.17) 



\u h (t) - u(t)\\ <Ch 2 £ h r a \\v\ 
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Proof. By Duhamel's principle 

5(t) = f F h (t - s)K h Q h d?u h (s)ds. 
Jo 

Then by appealing to the smoothing property of the operator F^ in Lemma 



4.2 and the inverse inequality (4.8), we get for X G ^h-, e > 0, and p = 0, 1 



(4.18) 



\\\F h (t)A h Q hX \\\ p < Ct^lWAhQhxlWp^+e 
= Ct^- 1 \\\Q hX \\\ p+e 

< Ct^- l h-%\Q hX \\\ p 

< Ct^h-'WQhxWp. 



Consequently, by Lemmas 



4.3 



3.1 



and H 1 - and ^-stability of the operator 



Ph from Lemma 2.2 we deduce for q = 0, 1 

\\V5(t)\\ <C^ +1 ~ £ [ (t-a)i a - x \\a?u h (a)\\ q d8 
Jo 

<Chi +1 ^ [ (t-8)%'*- 1 8- a da\\u h (p)\U 
Jo 

= Ch q+1 - £ r a ^-i)B 1 - a) \\P h v\ 

< c-h q+1 - e r a{l -^\v\ n . 



Now the desired estimate (4.14) follows by triangle inequality from the es- 
timate (3.23) and the above estimate by taking e = 1 and e = 1/lh for the 
cases q = 1 and 0, respectively. 

Next we derive an L2- error estimate. First, note that for X £ Xh we have 



\\F h (t)A h Q hX \\ < Ct^~ l \\\A h Q hX \\\-i < Ct^ L \\VQ hX \\. 



This estimate and Lemma 4.3 



give 



\W)\\ <Ch" +1 / (t-s^-^Uhis^ds 
Jo 

<Ch q+l f (t-8)?- 1 8- a d8\u h (0)\ q 

Jo 



<Ch q+l t-^B(^,l-a)\P h v\ q 
<Ch q+l t-1\v\ q , g = 0,l, 



which shows the desired estimate (4.15). 
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get 



Finally, if (4.16) holds, by applying (4.18) with p = and e G (0, \), we 

(t-s)i a - x \\Q h d?u h {s)\\ds 



\\5(t)\\<Ch- 

< Ch 2 - 

< Ch 2 - 
,1 



o 



(t 



S 2 



a-li 







(t 



*ds\u h (0)\ds 



<C ^ h 2-e t -a(l-l) l 



Then (4.17) follows immediately by choosing e = 1/lh- 



□ 



Remark 4.3. The condition (4.16) on the quadrature error operator Qy, 



is satisfied for symmetric meshes; see [2, Sections 5]. In case the condi- 



tion (4.16) does not hold, we were able to show only a suboptimal 0(h)- 



convergence rate for L2-norm of the error, which is reminiscent of the situ- 
ation in the classical parabolic equation (see, e.g. [2, Theorem 4.4]). 



Remark 4.4. As we mentioned before, assumption (4.16) is valid for sym 



metric meshes. In fact, in one dimension, the symmetry requirement can 



relaxed to almost symmetry [21 Section 6], and (4.17) can be proven as well. 

Remark 4.5. We note that we have used a globally quasi-uniform meshes, 
while the results in [2] are valid for meshes that satisfy the inverse inequality 
only locally. 

5. Special meshes 



Remark |3.3| (as well as Remark |4.1[ ) suggests that one can achieve a higher 
order convergence rate for V(it/i — u), if one can get an estimate of the 
error V(Rh,u — u) in some special norm. This could be achieved using the 
super convergence property of the gradient available for special meshes and 
solutions in H 3 (Q). Examples of special meshes exhibiting superconvergence 
property include triangulations in which every two adjacent triangles form a 
parallelogram [13J. To establish a super-convergent recovery of the gradient, 
Kfizek and Neittaanmaki in [J5] introduced an operator of the averaged 
(recovered, postprocessed) gradient Gh(Rhu) of the Ritz projection R^u of 
a function u (see \YA\ equation (2.2)]) with the following properties: 

(a) If u G H 3 (n) then, [E2 Theorem 4.2] 

(5.1) ||V« - G h (R h u)\\ < Ch 2 \\u\\ H3{n) . 

(b) For x G Xh the following bound is valid: 
(5-2) \\G h ( X )\\<C\\V X \\. 



The bound (5.2) follows immediately from |13| inequality (3.4)] established 
for a reference finite element by rescaling and using the fact that x £ Xh- 
We point out that one can get a higher order approximation of the Vw by 
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Gh(R-hu) due to the special post-processing procedure valid for sufficiently 
smooth solutions and special meshes. 

This result could be used to establish a higher convergence rate for the 
semi-discrete Galerkin method (and similarly for the lumped mass method) 
for smooth initial data. Specifically, we have the following result. 

Theorem 5.1. Let Th be strongly uniform triangulation ofVt, that is, every 
two adjacent triangles form a parallelogram. Then the following estimate is 
valid 



(5.3) 



\Vu(t)-G h {u h {t))\\<Ch 2 t- a l 2 \\v\ 



Proof. It follows from the fact that u satisfies equation ( |1.1[ ), i.e., dfu{t) = 
Au and from Theorem 2.1 (with £ = 1, p = 1 and q = 2) that 

Ht)\ 3 < Ct- a ' 2 \v\2. 

Then using the above super-convergent recovery operator Gh of the gradient 
with the properties (5.1), (5.2) and the estimate (3.19) for 9{t) = Rhu(t) — 
Uh(t), we get 

||V«(t) - G h (u h (t))\\ = \\Vu(t) - G h (R h u(t))\\ + \\G h (R h u(t) - u h {t))\\ 

<ch 2 \\u{t)\\ mn) + c\\ve{t)\\ 
< ch 2 r a ' 2 \\vh 

which shows the desired estimate. □ 



Remark 5.1. By repeating the proof of Theorem 5.1 and appealing to Re- 
mark 4-1, we can derive the following error estimate for the solution of the 
lumped mass Galerkin FEM 

\\Vu(t) - G h (u h (t))\\ < Ch 2 £ h t- a / 2 \\v\\ 2 , 4 = I ln/i|. 

Remark 5.2. Obviously any strongly regular triangulation is also symmetric 
at each internal vertex and therefore for such meshes we have as well optimal 
convergence in L2-norm for nonsmooth data; see (4.17). 



6. Numerical results 

In this section, we present some numerical results to verify the error es- 
timates. We consider the following one-dimensional problem on the unit 
interval (0, 1) 

d?u(x, t) - d 2 xx u{x, t) = 0, < x < 1 < t < T, 
(6.1) «(0, t) = u(l,t) = 0, < t < T, 

u(x,0) = v(x), < x < 1. 

We performed numerical tests on five different data: 



22 



BANGTI JIN, RAYTCHO LAZAROV, AND ZHI ZHOU 



(a) Smooth initial data: v(x) = — 4x 2 + 4x; in this case the initial data v 
is in H 2 (Q)nHQ (Q), and the exact solution u(x, t) can be represented 
by a rapidly converging Fourier series: 



u{x, t) = ^ V i_s (_ n Vr)(l - (-1)") 

n=l 

(b) Initial data in H 1 (intermediate smoothness): 
3.2) v(x) = 



sm nirx. 



xe [0,§], 

x, a; G (|, 1]. 



(c) Nonsmooth initial data: (1) u (x) = 1, (2) u(x) = x, and (3) v(x) = 
Xroii) the characteristic function of the interval (0, g)- Since this 
choice of v is not compatible with the homogeneous Dirichlet bound- 
ary data, obviously, in all three cases v £ Hq. However, in all these 
examples v G H s , < s < \. 

(d) We also consider initial data v that is a Dirac 5i (x)-function con- 

2 

centrated at x = |. Such weak data is not covered by our theory. 
However, it is interesting to see how the method performs for such 
highly nonsmooth initial data. We note that the choice of the Dirac 
delta function as initial data is common for certain parameter iden- 
tification problems in fractional diffusion problems [3]. 



(e) Variable coefficient case (cf. (3.24)): we take k(x) = 3 + sin(27rx) 



and initial condition v(x) = 1. This class of problems was discussed 
in Section [331 



The exact solution for each example from (a) to (d) can be expressed by an 
infinite series involving the Mittag-Leffler function E a \(z). To accurately 
evaluate the Mittag-Leffler functions, we employ the algorithm developed 
in [23], which is based on three different approximations of the function, 
i.e., Taylor series, integral representations and exponential asymptotics, in 
different regions of the domain. 

We divide the unit interval (0, 1) into N + 1 equally spaced subintervals, 
with a mesh size h = 1/(N + 1). The space consists of continuous piece- 
wise linear functions on the partition. In the case of a constant coefficient 
k(x) (cf. Section |3.4|) we can represent the exact solution to the semidis- 



crete problem by (3.2) for the standard semidiscrete Galerkin method and 
by (4.5) for the lumped mass method using the eigenpairs (Xj,ipj(x)) and 
(\j,(fij(x)) of the respective one-dimensional discrete Laplacian — and 
-Aft, i.e., 

(-A h ^,v) = X^^,v) and (-A h $>J, v) h = Aj($, v) h Vv G X h . 
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Here (w,v) and (w,v)h refer to the standard L2-inner product and the ap- 



proximate L2-inner product (4.2) on the space Xh, respectively. Then 



X j = X j = p sin2 2 (N + 1) and = ^i( Xfc ) = V2sm(jirx k ), j = 1,2,. 

for Xk being a mesh point and linear over the finite elements. These will be 
used in computing the finite element approximations by the Galerkin and 
lumped mass methods. 

We also have used a direct numerical solution technique by first discretiz- 
ing the time interval, t n = tit, n = 0, 1, . . . , with r being the time step size, 
and then using a weighted finite difference approximation of the fractional 
derivative dfu{x,t n ) developed in [16] : 



1 n_1 f 



" +1 duM (t n -s)-«ds 



3=0 
n-1 



ds 



r( 



1 — a) ' r ./*, 



3=0 

1 , "(a:, tn-j) - U(x, tn-j-l) 



r(2-«; 



E 6 

3=0 



where the weights bj = (j + l) 1_a — j l ~ a , j = 0, 1, . . . , n — 1. It has been 
shown that if the solution u(x, t) is sufficiently smooth and the time step r is 
a constant, then local truncation error of the approximation is bounded by 
Cr 2 ~ a for some C depending only on u [TBI equation (3.3)]. When needed, 
we have used this approximation on very fine meshes in both space and 
time to compute a reference solution. Unless otherwise specified, we have 
set r = 1.0 x 10~ 6 , so that the error incurred by temporal discretization 
can be ignored. Whenever possible, we have compared the accuracy of this 
reference solution with the exact representation. Our experiments show 
that with a very small time step size, these two produce the same numerical 
results. 

For each example, we measure the accuracy of the approximation Uh(t) 
by the normalized error \\u(t) — itft(£)[|/|M| and \\d x (u(t) — Uh(t))\\/\\v\\. The 
normalization enables us to observe the behavior of the error with respect 
to time in case of nonsmooth initial data. We shall present only numerical 
results for the lumped mass FEM, since that for the Galerkin FEM is almost 
identical. 



6.1. Numerical experiments for the smooth initial data: exam- 
ple (a). In Table [l] we report the numerical results for t = 1 and a = 
0.1, 0.5, 0.95. In Figure [TJ we show plots of the results from Table [JJ in 
a log-log scale. We see that the slopes of the error curves are 2 and 1, re- 
spectively, for L/2- and ff 1 -norm of the error. In the last column of Table 
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[I] we also present the error of the recovered gradient Gh(uh)- Since in one- 
dimension the mid-point of each interval has the desired superconvergence 
property, the recovered gradient in the case is very simple, just sampled 
at these points |26l Theorem 1.5.1]. It is clear that the recovered gradient 
Gh(uh) exhibits a 0(h 2 ) convergence rate, concurring with the estimates in 
Theorem |5.1| and Remark |5.1[ It is worth noting that the smaller is the a 
value, the larger is the error (in either the L2- or i? 1 -norm). This is at- 
tributed to the property of the Mittag-Leffler function E a ^(— Xt a ), which, 
asymptotically, decays faster as a approaches unity, cf . Lemma |2.1| and the 



representation (2.3). 



Table 1. Numerical results for smooth data, example (a): 
t = 1 and h = 2~ k . 





L2-error 


ZZ^-error 


Gh(uh)-evTOT 


k 


a = 0.1 


a = 0.5 


a = 0.95 


a = 0.1 


a = 0.5 


a = 0.95 


a = 0.5 


3 


5.23e-4 


3.37e-4 


4.84e-5 


2.65e-2 


1.74e-2 


2.04e-3 


3.20e-3 


4 


1.29e-4 


8.31e-5 


1.21e-5 


1.33e-2 


8.77e-3 


1.02e-3 


8.07e-4 


5 


3.21e-5 


2.07e-5 


3.05e-6 


6.69e-3 


4.39e-3 


5.11e-4 


2.03e-4 


6 


8.01e-6 


5.17e-6 


7.93e-7 


3.34e-3 


2.19e-3 


2.55e-4 


5.17e-5 


7 


2.00e-6 


1.30e-6 


2.32e-7 


1.67e-3 


1.10e-3 


1.28e-4 


1.39e-5 






~ L 2 


a=0.1 




— "H 


,a=0.1 




— L 2 


a-0.5 




— "H 


,a=0.5 






a=0.95 




-,- H 


,a=0.95 






(u h ),a=0.5 



10 

error 



Figure 1. Numerical results for smooth initial data, exam- 
ple (a) with a = 0.1,0.5,0.95 at t = 1. 



6.2. Numerical experiments for the intermediate case of smooth- 
ness of the data, example (b). In this example the initial data v(x) is 
in Hq(£1) n i72~ e (ri) with e > 0, and thus it represents an intermediate 
case. All the numerical results reported in Table [2] were evaluated at t = 1 
for a = 0.5. The slopes of the error curves in a log- log scale are 2 and 1 
respectively for L2 and H 1 norm of the errors, which is in agreement with 
the theory for the intermediate case; ratio in the last column of Table [2j 
refers to the ratio between the errors as the mesh size h halves. 
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Table 2. Numerical results for the intermediate case (b) 
with a = 0.5 at t = 1. 



h 


1/8 1/16 1/32 1/64 1/128 


ratio 


L2-error 
.f/^-error 


8.08e-4 2.00e-4 5.00e-5 1.26e-5 3.24e-6 
1.80e-2 8.84e-3 4.39e-3 2.19e-3 1.10e-3 


« 3.97 
« 2.00 



6.3. Numerical experiments for nonsmooth initial data: example 
(c). In Tables [3] and [4] we present the computational results for problem (c), 
cases (1) and (2). For nonsmooth initial data, we are particularly interested 
in errors for t close to zero, and thus we also present the error at t = 0.005 
and t = 0.01. In Figure [2] we plot the results shown in Tables [3] and |4j i.e., 
for problem (c), cases (1) and (2). These numerical results fully confirm the 
theoretically predicted rates for the nonsmooth initial data. 




10" 









^L 2 ,t=0.005 




1 s ?s : 
1 ^ / // 




— H',t=0.005 
^_L 2 ,t-0.01 










-— H',t=0.01 


JlC^ / I 
















H',W 


I I I I I 







io" b io" 5 10" 4 io" 3 10" Z 10" 



error 

(b) Error plots for Example (c) (2) 



Figure 2. Numerical results for nonsmooth initial data with 
a = 0.5 at t = 0.005, 0.01, 1.0. 

Now we consider the third example of nonsmooth case, the characteristic 
function of the interval (0,0.5), namely, v(x) = X\o ii- Note that if we use 
the interpolation of v as the initial data for the semidiscrete problem, the L^- 
error has only a suboptimal first-order convergence. However, if we choose 
L>2 projection as is discussed in previous sections, then the results agree well 
with our estimates; see Table [5] We also discretize this example by the 
Galerkin method, and the results are presented in Table [6] A comparison of 
Tables [5] and [6] clearly indicates that the Galerkin method and the lumped 
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Table 3. Nonsmooth initial data, example (c) (1): a = 0.5 
at t = 0.005,0.01,1. 





h 


1/8 1/16 1/32 1/64 1/128 


ratio 


t = 0.005 


L2-norm 


1.06e-2 2.65e-3 6.63e-4 1.65e-4 4.02e-5 


ps 4.05 


ii^-norm 


2.08e-l 1.04e-l 5.22e-2 2.61e-2 1.30e-2 


ps 2.00 


t = 0.01 


L2-norm 


7.94e-3 1.99e-3 4.93e-4 1.19e-4 2.59e-5 


» 4.08 


ii^-norm 


1.63e-l 8.16e-2 4.08e-2 2.04e-2 1.02e-2 


ps 2.00 


t = 1 


L2-norm 


8.07e-4 2.02e-4 5.03e-5 1.25e-5 3.05e-6 


ps 4.02 


ii^-norm 


2.02e-2 1.01e-2 5.04e-3 2.52e-3 1.26e-3 


ps 2.00 


Table 4. Nonsmooth initial data, example (c)(2): a = 0.5 
at t = 0.005,0.01,1. 


Time 


h 


1/8 1/16 1/32 1/64 1/128 


ratio 


t = 0.005 


L2-norm 


1.08e-2 2.71e-3 6.79e-4 1.69e-4 4.13e-5 


» 4.03 


if 1 -norm 


2.28e-l 1.14e-2 5.71e-2 2.86e-2 1.43e-2 


« 2.00 


t = 0.01 


L2-norm 


7.98e-3 2.00e-3 4.99e-4 1.23e-4 2.91e-5 


» 4.02 


ii^-norm 


1.73e-l 8.67e-2 4.34e-2 2.17e-2 1.08e-2 


» 2.00 


t = 1 


L2-norm 


8.05e-4 2.01e-4 5.03e-5 1.25e-5 3.07e-6 


» 4.01 


ii^-norm 


2.02e-2 1.01e-2 5.07e-3 2.53e-3 1.27e-3 


« 2.00 



mass method yield almost identical results for this example. Although not 
presented, we note that similar observations hold for other examples as well. 
Hence, we have focused our presentation on the lumped mass method. 

Table 5. Nonsmooth initial data, example (c)(3): a = 0.5 
at t = 0.005,0.01,1 



time 


h 


1/8 1/16 1/32 1/64 1/128 


ratio 


t = 0.005 


L2-norm 


8.54e-3 2.16e-3 5.45e-4 1.31e-4 3.17e-5 


ps 4.06 


i/ 1 -norm 


2.18c-l 1.08e-l 5.38e-2 2.68e-2 1.34e-2 


ps 2.00 


t = 0.01 


L2-norm 


6.54e-3 1.64e-3 4.14e-4 1.06e-4 2.90e-5 


ps 3.96 


i/ 1 -norm 


1.63e-l 8.04e-2 4.00e-2 2.00e-2 9.96e-3 


ps 2.00 


t = 1 


L2-norm 


8.10e-4 2.03e-4 5.07e-5 1.27e-5 3.23e-6 


ps 3.99 


i/ 1 -norm 


1.82e-2 9.02e-3 4.46e-3 2.22e-3 l.lle-3 


ps 2.01 



6.4. Numerical experiments for initial data a Dirac (5-function. We 

note that this case is not covered by our theory. Formally, the orthogonal L2- 
projection is not well defined for such functions. However, we can look at 
(v, x) for x € -X/i C Hq(Q) as a duality pairing between the spaces i? _1 (0) 
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Table 6. Nonsmooth initial data, example (c)(3) with a = 
0.5 at t = 0.005, 0.01, 1 by the Galerkin method. 



time 


h 


1/8 1/16 1/32 1/64 1/128 


ratio 


t = 0.005 


L2-norm 


8.60e-3 2.14e-3 5.30e-4 1.28e-4 2.85e-5 


ps 4.11 


iT -norm 


1.78e-l 9.78e-2 5.11e-2 2.61e-2 1.32e-2 


ps 1.96 


t = 0.01 


L2-norm 


6.56e-3 1.64e-3 4.06e-4 9.94e-5 2.29e-5 


ps 4.11 


iJ -norm 


1.34e-l 7.34e-2 3.82e-2 1.95e-2 9.85e-3 


ps 1.95 


t = 1 


L2-norm 


8.07e-4 2.02e-4 5.04e-5 1.25e-5 3.09e-6 


ps 4.03 


iJ -norm 


1.54e-2 8.30e-3 4.29e-3 2.18e-3 1.10e-3 


ps 1.96 



and Hq(Q) and therefore (5i,x) = x{\)- If 

mesh point, say 

xl, then we can define Ph&i appropriately with its finite element expansion 

2 

given by the Lth column of the inverse of the mass matrix. This was the 
initial data for the semidiscrete problem that we used in our computations. 
In Table [7] we show the Li- and // 1 -norm of the error for this case. It is 
noted that the i7 1 -norm of the error converges as 0(h%), while the error in 
the L2-norm converges as O(h^); see the last column of Table [Tj It is quite 
remarkable that we can practically have good convergence rates in L2- and 
i? 1 -norm for such very weak solutions. A theoretical justification of these 
rates is a subject of our current work. 

Table 7. Lumped mass FEM with initial data a Dirac 5- 
function, a = 0.5, t = 0.005, 0.01, 1. 



time 


h 


1/8 1/16 1/32 1/64 1/128 


ratio 


t = 0.005 


L2-norm 


7.24e-2 2.66e-2 9.54e-3 3.40e-3 1.21e-3 


ps 2.79 


i? 1 -norm 


1.51e0 1.07e0 7.60e-l 5.40e-l 3.81e-l 


ps 1.41 


t = 0.01 


L2-norm 


5.20e-2 1.89e-2 6.77e-3 2.40e-3 8.54e-4 


ps 2.79 


ff i -norm 


1.07e0 7.59e-l 5.37e-l 3.80e-l 2.70e-l 


ps 1.41 


t = 1 


L2-norm 


5.47e-3 1.93e-3 6.84e-4 2.42e-4 8.56e-5 


ps 2.79 


ff 1 -norm 


1.07e-l 7.58e-2 5.37e-2 3.80e-2 2.70e-2 


ps 1.41 



6.5. Numerical experiments for variable coefficients, example (e). 

Although we do not have an explicit representation of the exact solution, 
we compare the numerical solution with the approximate solution obtained 
on very fine meshes, namely, with mesh-size h = 1/512 and time-step size 
r = 1.0 x 10 -5 . Normalized L2- and i? 1 -norms of the error are reported 
in Table [8] for t = 0.01 for a = 0.5. The results confirm the theoretically 
predicted rate. 

In summary, the convergence rates observed for all the numerical ex- 
periments are in excellent agreement with the theoretical findings for both 
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Table 8 . Numerical results for variable coefficients and non- 
smooth initial data with a = 0.5 at t = 0.01. 



h 


1/8 1/16 1/32 1/64 1/128 


ratio 


L2-error 
f/^-error 


3.24e-3 8.21e-4 2.05e-4 5.09e-5 1.23e-6 
7.15e-2 3.60e-2 1.80e-2 8.94e-3 4.36e-3 


ps 4.02 
« 2.01 



smooth and nonsmooth initial data, including also the case of the recovered 
gradient Gh(uh) discussed in Section [5j see also Remark 5.1 
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